"""
数据拟合与插值算法

问题描述：
    数据拟合与插值是数据处理中的基础技术，用于从有限的数据点中推断连续函数。
    拟合是寻找一个简单函数来近似描述数据的总体趋势，允许一定误差；
    插值是构造一个精确通过所有数据点的函数。

算法原理：
    1. 多项式拟合：使用多项式函数逼近数据
    2. 样条插值：分段多项式插值，保证光滑性
    3. 最小二乘拟合：最小化残差平方和
    4. 核函数拟合：使用径向基函数等核函数
    5. 贝塞尔曲线：参数化曲线拟合

核心概念：
    - 拟合误差：数据点与拟合函数的偏差
    - 光滑性：函数的连续可导性
    - 过拟合：模型过于复杂，泛化能力差
    - 交叉验证：评估拟合效果的方法

应用场景：
    数据拟合与插值广泛应用于：
    - 实验数据分析
    - 数值计算
    - 图像处理
    - 信号处理
    - 工程设计

作者：斯黄
日期：2025年3月4日
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from scipy.optimize import curve_fit
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
import warnings
warnings.filterwarnings('ignore')

class DataFitting:
    def __init__(self):
        self.fitted_functions = {}
        
    def polynomial_fitting(self, x, y, degree=3):
        """多项式拟合"""
        coefficients = np.polyfit(x, y, degree)
        poly_func = np.poly1d(coefficients)
        
        # 计算拟合优度
        y_pred = poly_func(x)
        r_squared = 1 - np.sum((y - y_pred)**2) / np.sum((y - np.mean(y))**2)
        
        self.fitted_functions['polynomial'] = {
            'function': poly_func,
            'coefficients': coefficients,
            'degree': degree,
            'r_squared': r_squared
        }
        
        return poly_func, r_squared
    
    def cubic_spline_interpolation(self, x, y, kind='cubic'):
        """样条插值"""
        spline_func = interpolate.interp1d(x, y, kind=kind, fill_value='extrapolate')
        
        self.fitted_functions['spline'] = {
            'function': spline_func,
            'kind': kind
        }
        
        return spline_func
    
    def gaussian_fitting(self, x, y):
        """高斯函数拟合"""
        def gaussian(x, a, mu, sigma):
            return a * np.exp(-(x - mu)**2 / (2 * sigma**2))
        
        # 初始参数估计
        a_init = np.max(y)
        mu_init = x[np.argmax(y)]
        sigma_init = (np.max(x) - np.min(x)) / 4
        
        try:
            popt, pcov = curve_fit(gaussian, x, y, p0=[a_init, mu_init, sigma_init])
            
            def fitted_gaussian(x_new):
                return gaussian(x_new, *popt)
            
            # 计算拟合优度
            y_pred = fitted_gaussian(x)
            r_squared = 1 - np.sum((y - y_pred)**2) / np.sum((y - np.mean(y))**2)
            
            self.fitted_functions['gaussian'] = {
                'function': fitted_gaussian,
                'parameters': popt,
                'covariance': pcov,
                'r_squared': r_squared
            }
            
            return fitted_gaussian, popt, r_squared
        except:
            print("高斯拟合失败")
            return None, None, None
    
    def exponential_fitting(self, x, y):
        """指数函数拟合"""
        def exponential(x, a, b, c):
            return a * np.exp(b * x) + c
        
        # 初始参数估计
        a_init = (y[-1] - y[0]) / (np.exp(x[-1]) - np.exp(x[0]))
        b_init = 1.0
        c_init = y[0]
        
        try:
            popt, pcov = curve_fit(exponential, x, y, p0=[a_init, b_init, c_init])
            
            def fitted_exponential(x_new):
                return exponential(x_new, *popt)
            
            # 计算拟合优度
            y_pred = fitted_exponential(x)
            r_squared = 1 - np.sum((y - y_pred)**2) / np.sum((y - np.mean(y))**2)
            
            self.fitted_functions['exponential'] = {
                'function': fitted_exponential,
                'parameters': popt,
                'covariance': pcov,
                'r_squared': r_squared
            }
            
            return fitted_exponential, popt, r_squared
        except:
            print("指数拟合失败")
            return None, None, None
    
    def fourier_fitting(self, x, y, n_harmonics=3):
        """傅里叶级数拟合"""
        def fourier_series(x, *coeffs):
            result = coeffs[0]  # 常数项
            for i in range(1, len(coeffs)//2 + 1):
                a_n = coeffs[2*i-1]
                b_n = coeffs[2*i]
                result += a_n * np.cos(2*np.pi*i*x/(x[-1]-x[0])) + b_n * np.sin(2*np.pi*i*x/(x[-1]-x[0]))
            return result
        
        # 初始参数（常数项 + n_harmonics对正弦余弦系数）
        p0 = [np.mean(y)] + [0.1] * (2 * n_harmonics)
        
        try:
            popt, pcov = curve_fit(fourier_series, x, y, p0=p0)
            
            def fitted_fourier(x_new):
                return fourier_series(x_new, *popt)
            
            # 计算拟合优度
            y_pred = fitted_fourier(x)
            r_squared = 1 - np.sum((y - y_pred)**2) / np.sum((y - np.mean(y))**2)
            
            self.fitted_functions['fourier'] = {
                'function': fitted_fourier,
                'parameters': popt,
                'n_harmonics': n_harmonics,
                'r_squared': r_squared
            }
            
            return fitted_fourier, popt, r_squared
        except:
            print("傅里叶拟合失败")
            return None, None, None
    
    def rbf_interpolation(self, x, y, function='multiquadric'):
        """径向基函数插值"""
        # 创建二维数据点（假设是一维数据）
        points = np.column_stack([x, np.zeros(len(x))])
        values = y
        
        rbf = interpolate.Rbf(x, y, function=function)
        
        self.fitted_functions['rbf'] = {
            'function': rbf,
            'rbf_function': function
        }
        
        return rbf

class InterpolationMethods:
    @staticmethod
    def lagrange_interpolation(x_data, y_data, x_new):
        """拉格朗日插值"""
        n = len(x_data)
        y_new = np.zeros_like(x_new)
        
        for i in range(n):
            # 计算拉格朗日基函数
            L_i = np.ones_like(x_new)
            for j in range(n):
                if i != j:
                    L_i *= (x_new - x_data[j]) / (x_data[i] - x_data[j])
            y_new += y_data[i] * L_i
        
        return y_new
    
    @staticmethod
    def newton_interpolation(x_data, y_data, x_new):
        """牛顿插值"""
        n = len(x_data)
        
        # 计算差商表
        divided_diff = np.zeros([n, n])
        divided_diff[:, 0] = y_data
        
        for j in range(1, n):
            for i in range(n - j):
                divided_diff[i, j] = (divided_diff[i+1, j-1] - divided_diff[i, j-1]) / (x_data[i+j] - x_data[i])
        
        # 计算插值结果
        y_new = np.zeros_like(x_new)
        for i in range(n):
            term = divided_diff[0, i]
            for j in range(i):
                term *= (x_new - x_data[j])
            y_new += term
        
        return y_new

def demonstrate_fitting():
    """演示拟合方法"""
    print("=== 数据拟合演示 ===")
    
    # 生成测试数据
    np.random.seed(42)
    x = np.linspace(0, 10, 20)
    y = 2 * np.sin(x) + 0.5 * x + np.random.normal(0, 0.3, len(x))
    
    fitter = DataFitting()
    
    # 多项式拟合
    poly_func, r2 = fitter.polynomial_fitting(x, y, degree=4)
    print(f"多项式拟合 R² = {r2:.4f}")
    
    # 可视化
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    x_dense = np.linspace(0, 10, 200)
    plt.scatter(x, y, color='red', s=50, label='原始数据')
    plt.plot(x_dense, poly_func(x_dense), 'b-', linewidth=2, label='多项式拟合')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('多项式拟合')
    plt.legend()
    plt.grid(True)
    
    # 样条插值
    spline_func = fitter.cubic_spline_interpolation(x, y)
    plt.subplot(1, 2, 2)
    plt.scatter(x, y, color='red', s=50, label='原始数据')
    plt.plot(x_dense, spline_func(x_dense), 'g-', linewidth=2, label='三次样条')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('样条插值')
    plt.legend()
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    demonstrate_fitting() 